# -*- coding: utf-8 -*-
"""
Created on Wed Jan 27 13:22:09 2021

@author: Lenovo
"""

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
# code for the MCMC algorithm

from scipy.stats import norm

# define a sum of squares function that we will use as fitness
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import leastsq

#Load the data
didata1 = np.loadtxt("amrdata.txt")
time = didata1[:,0]
pop = didata1[:,1]

# defines the model equations
def diauxic_ode(x,t,params):
    r1,r2,k = params  
    y,S1,S2 = x
    derivs = [r1*S1*y+(k/(k+S1))*r2*S2*y, -r1*S1*y, -(k/(k+S1))*r2*S2*y]
    return derivs

# runs a simulation and returns the population size
def diauxic_run(pars,t):
    r1,r2,k,y0,S10,S20 = pars
    ode_params=[r1,r2,k]
    ode_starts=[y0,S10,S20]
    out = odeint(diauxic_ode, ode_starts, t, args=(ode_params,))
    return out[:,0] # we only return the population size - we don't worry about the substrates as they are not measured


#log likelihood fuction for the diauxic equation
def diauxic_loglik(pars,t,data,sig):
    out = diauxic_run(pars,t)
    return sum(norm.logpdf(out,data,sig))



# How many replicate simulations will you need? More than 1.
reps = 1

# parameters for the MCMC
npars = 6

# set standard deviations to be 10% of the population values
sigma = pop/10

# output matrix
par_out = np.ones(shape=(reps,npars))

# 通过对初值的一些好的猜测来完成这段代码
par_out[0,] = [1.0,1.0,1.0,1.0,1.0,1.0]

# plot a run the first set of parameters
plt.plot(time,pop,'h')
plt.show()

plt.plot(time,diauxic_run(par_out[0,:],time))
plt.show()

# acceptance 
accept = np.zeros(shape=(reps,npars))

# 建议标准差。根据需要调整这些
propsigma = [0.1,0.1,0.1,0.1,0.1,0.1]

for i in range(1,reps):
    # make a copy of previous parameters
    par_out[i,] = par_out[i-1,]
    for j in range(npars):
        proposed = np.copy(par_out[i,:]) # we need to make a copy so that rejected moves don't affect the original matrix
        # use a log normal proposal
        proposed[j] = np.exp(np.log(proposed[j]) + np.random.normal(0,propsigma[j]))
        alpha = np.exp(diauxic_loglik(proposed,time,pop,sigma)-diauxic_loglik(par_out[i,],time,pop,sigma))
        u = np.random.rand()
        if (u < alpha):
            par_out[i,j] = proposed[j]
            accept[i,j] = 1


# plot graph of last sample - this is just to check it is working OK
plt.plot(time,diauxic_run(par_out[-1,:],time))
plt.show()